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Abstract. - A fitting scheme is proposed to obtain effective potentials from Car-Parrinello 
molecular dynamics (CPMD) simulations. It is used to parameterize a new pair potential for 
silica. MD simulations with this new potential are done to determine structural and dynamic 
properties and to compare these properties to those obtained from CPMD and a MD simulation 
using the so-called BKS potential. The new potential reproduces accurately the liquid structure 
generated by the CPMD trajectories, the experimental activation energies for the self-diffusion 
constants and the experimental density of amorphous silica. Also lattice parameters and elastic 
■ constants of a— quartz are well-reproduced, showing the transferability of the new potential. 
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Introduction. — One of the central issues of atom- 
istic computer simulations is to develop effective poten- 
tials [1] that model the interactions between ions by some 
functional form without taking into account explicitly the 
electronic degrees of freedom. The first step to obtain such 
a potential for a given system is to decide its functional 
form. In terms of computational efficiency, an optimal 
ansatz is provided by a pair potential model. However, 
for a realistic modelling it might be often necessary to 
include many-body terms to describe, e.g., three-body in- 
teractions or polarizability effects. In a second step, the 
free parameters of the potential function have to be fitted 
to experimental data and/or to ab initio calculations [2,3]. 
A relatively new approach is to use ab initio molecular dy- 
namics simulations such as Car-Parrinello molecular dy- 
namics (CPMD) [4] for the parametrization of effective po- 
tentials. The latter simulation technique has been proven 
very successful to realistically describe a large number of 
atomistic systems, among them network forming systems 
such as silica [5-10]. Nowadays, systems of 100 to 200 
particles can be simulated on a time scale of several tens 



of picoseconds by CPMD. 

In this Letter, we aim at developing a fitting scheme 
for the parametrization of effective potentials from CPMD 
trajectories and we address this issue to the case of liquid 
silica. A widely-used method for the parametrization of 
potentials is the so-called force matching procedure [11] 
where one tries to match the CPMD forces on the atoms 
with those resulting from the classical potential. How- 
ever, we have found that, in the case of silica, the force 
matching method fails (at least if a pair potential model 
is considered) [12]: Although the CPMD forces are well 
reproduced by the resulting effective (pair) potential, the 
structure is very different from that obtained by CPMD. 
The main aim of this work is to develop a different fitting 
scheme where one tries to match as closely as possible the 
atomistic structure, as obtained from the classical MD us- 
ing the effective potential, with that of the CPMD. 

Since silica is of great importance for technological ap- 
plications, geosciences and the theory of the glass tran- 
sition, the development of effective potentials for silica 
has a long history [2,3,13-18]. Surprisingly, various stud- 
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ies have shown that a simple pair potential, the so-called 
BKS potential [2] , is able to give a quite accurate descrip- 
tion of amorphous silica [8,19-22]. The parametrization 
of this potential was based on Hartree-Fock calculations 
of a single Si04 tetrahedron (saturated by four hydrogen 
atoms), considering also experimental clastic constants of 
a— quartz in the fitting procedure. 

Unsurprisingly, the BKS model has also several defi- 
ciencies. For instance, the equation of state, as far as it 
is known experimentally, is not reproduced well [3]. One 
could argue that one has to include many body effects to 
overcome the deficiencies of the BKS potential. However, 
as shown in Rcf. [22], more complicated potentials that 
include polarizability effects or fluctuating charges do not 
provide necessarily a more realistic description of silica. 

Here, we demonstrate that one can improve the accu- 
racy of a pair potential for amorphous silica. The idea of 
our approach is to match the partial pair correlation func- 
tions, as obtained from the effective potential, with those 
obtained from a CPMD simulation. As demonstrated be- 
low, the new potential yields an accurate description of 
amorphous silica with respect to density, structure and 
diffusion dynamics. Also properties of a— quartz are well 
reproduced by the new potential. The fitting scheme pro- 
posed in this work can be also used for other systems and it 
is not restricted to the parametrization of pair potentials. 

Simulation details and fitting methodology. 

Ab initio simulations. CPMD simulations [23] were 
done for a system of 38 SiC>2 units (114 atoms) in a cu- 
bic box of size L = 11.982 A, corresponding to the ex- 
perimental density of 2.2g/cm 3 [24]. Periodic boundary 
conditions are used in all three spatial directions. The 
electronic structure was described by the local density ap- 
proximation (LDA) in the framework of density functional 
theory [25]. Norm-conserving pseudopotcntials were used 
for silicon [26] and oxygen [27], and the plane- wave basis 
sets were expanded up to an energy cut-off of 70 Ry. The 
choice of pseudopotcntials, exchange and correlation func- 
tionals and the plane- wave cutoff are justified by previous 
CPMD studies carried out on amorphous SiC>2 [8-10]. 

A time step of 3 a.u. (0.0725 fs) and a fictitious electronic 
mass of 600 a.u. were used to integrate the equations of 
motion. Temperature is kept constant by applying Nose- 
Hoover chains for each ionic degree of freedom [28,29] as 
well as for the electronic degrees of freedom to counter- 
balance the energy flow from ions to electrons [30]. The 
parameters used for the Nose-Hoover chains can be found 
in previous publications [8,9]. 

The starting configuration for the CPMD was generated 
by MD simulations with the BKS potential [2], following 
the same methodology as in Refs. [8,9]. Then, a CPMD 
run over about 20 ps at T = 3600 K was performed from 
which we used the last 16.5 ps for the analysis. 

Fitting procedure. The next step was to obtain an ef- 
fective potential from the CPMD data. To this end, we 
aimed at parameterizing a given potential model such that 



the partial pair correlation functions g a p{r) (u(3 £ {Si, O}) 
[31], as calculated from the CPMD trajectories, were re- 
produced by the MD simulations using the new potential. 

As an ansatz for the effective potential we used the same 
functional form as the one used for the BKS potential, 

u a p{r) = — \- A a pcxjp{-B af} r) — , (1) 

where r is the distance between an ion of type a and 
an ion of type f3 {a, (3 = Si, O) and e is the elementary 
charge. We consider the parameters appearing in Eq. (1) 
as elements of the vector £ = k = 1, 2, . . . 10) = 
(qsi,A a {3,B a i3,C a f3,a,f3 £ {Si, O}). The effective charge 
for the oxygen atoms does not appear in the £ vector since, 
due to the charge neutrality requirement, it is directly re- 
lated to the silicon charge via go = — 9si/2- Forces and en- 
ergies that correspond to the long-ranged Coulomb terms 
in (1) have been computed by Ewald sums [32], those for 
the short-range part in (1) have been truncated and shifted 
at 6.5 A. In addition a smoothening function was used for 
the short-range part to avoid a drift of the total energy in 
microcanonical MD runs. The details of this function can 
be found in a recent MD study with the BKS model [33]. 

The cost function x 2 that measures the difference be- 
tween the CPMD data and the fitting model is based on 
the pair correlation functions weighted by the distance r, 

,L/2 

X 2 = / \ 2 (r)dr, (2) 
Jo 

where A 2 (r) is defined as follows: 

A 2 M= £ [rs^M-rj^rifl] 8 . (3) 

as,/3=Si,0 

Here, the superscripts CPMD and CHIK indicate respec- 
tively the use of CPMD and MD with the new potential, 
that we call CHIK potential in the following. The integral 
(2) was evaluated by numerical integration. Note that a 
cost function \ 2 m which the function rg a p(r) in Eq. (3) 
is replaced by gapif) or r 2 g a p(r) resulted in potentials 
that were not accurate [12]. The reason for this is that 
rg a p (r) makes a good compromise between weighting the 
structural information at short distances, such as is done 
by the cost function with g a /3(r), and long distances, as is 
done by the cost function with r 2 g a p(r). 

The optimal vector £ of the potential parameters was 
determined by an iterative Levenberg-Marquardt algo- 
rithm [34]. After each iteration step, £ was updated 
leading to a new guess for the potential. This poten- 
tial was then used in MD simulations yielding a new set 
of pair correlation functions ga^ lK ( r l £)■ Recall that the 
Levenberg-Marquardt algorithm requires the calculation 
of the derivatives of the g a p with respect to the potential 
parameters. These derivatives were computed numerically 
by finite differences: 

dg al 3(r,£) = ^ g a p(r, £ fc + e k ) - g a p{r, g fc - e k ) 
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Parameter 


CHIK 


<ZSi 


(C) 


1.910418 


Ado 


(eV) 
(A" 1 ) 


659.595398 


-Boo 


2.590066 


Coo 


(eV.A 6 ) 


26.836679 


A Si0 


(eV) 
(A" 1 ) 


27029.419922 


Bsio 


5.158606 


Csio 


(eV.A 6 ) 


148.099091 


A SiSi 


(eV) 


3150.462646 


B SiSi 


(A- 1 ) 


2.851451 


Csisi 


(eV.A 6 ) 


626.751953 



1.5 



Table 1: Parameters of the CHIK potential. 
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Here, small perturbations (±£fc) of a given parameter 
are considered, keeping all other parameters constant. 

MD simulations were performed to compute the func- 
tions 9a™ K ( r i an d the derivatives (4). In these sim- 
ulations, we considered a system of 1152 atoms in a 
cubic box of size 25.9 A (corresponding to the density 
p = 2.2 g/cm 3 ). The simulations were done in the NVT 
ensemble at the temperature T = 3600 K. Temperature 
was kept constant by coupling the system periodically to 
a stochastic heat bath. The equations of motion were inte- 
grated by the velocity form of the Verlet algorithm, using 
a time step of 1.0 fs. For an accurate calculation of the 
functions 5 , q™ k (t';0; runs over 10 ps were required. 

Following the methodology described above, conver- 
gence was obtained after 44 iterations yielding the poten- 
tial parameters 1 that are listed in Tab. 1. We note that 
these parameters are of the same order of magnitude as 
the parameters of the BKS potential [2] . Thus, our model 
can be seen as a modification of the BKS potential which 
yields already a good description of (amorphous) silica. 

MD simulations with the CHIK potential. As a third 
step extensive MD simulations were performed with the 
CHIK potential to investigate the temperature depen- 
dence of structural and dynamic properties of the new 
model. We used systems of 1152 particles in a cubic 
box of size L = 25.9 A corresponding to the density 
p = 2.2 g/cm 3 . Equilibration runs were done at con- 
stant temperature (see above), using now a time step 
of 1.6 fs. After the equilibration, microcanonical produc- 
tion runs were done at the temperatures 5200 K, 4300 K, 
4000 K, 3580 K, 3250 K, 3000 K, 2900 K, 2750 K, 2580 K, 
and 2440 K. At T = 2440 K, each the equilibration and 
production runs lasted over 32 million time step corre- 
sponding to a real time of 52.3 ns. At each temperature, 



*As the Born-Mayer ansatz diverges for very short distances, the 
following repulsive interaction has been added to the potential ansatz 
(1): 



with Doo = 113.0 eV. A 24 , D Si0 
3423200.0 eV. A 24 . 



29.0 eV.A 24 , Dsisi 



Fig. 1: Temperature dependence of the pressure for the BKS 
model (from Ref. [20]) and the CHIK model, as indicated. 



results were averaged over 8 independent runs. In order to 
estimate the pressure, additional runs were done at 2230 K 
and 2000 K. Here, the system was annealed for 15 million 
time steps. Then the pressure was measured in subsequent 
runs over 5 million time steps. 

For comparison, also results from a previous MD sim- 
ulation with the BKS model are shown in the following. 
The details of this simulation can be found in Ref. [20] . 

Results. 

Static properties. One of the peculiar properties of 
amorphous Si02 is the occurrence of a density maximum 
at 1820 K [35]. Moreover, in a broad temperature range 
from room temperature to temperatures above 2000 K the 
variations of the density around an average value of about 
2.2 g/cm 3 are smaller than 1% [24]. These features can 
also be seen in the temperature dependence of the pressure 
since the pressure shows a local minimum at the temper- 
ature where the density maximum occurs. Both the BKS 
and the CHIK model show indeed a minimum pressure at 
about 4800 K and 2300 K, respectively, see Fig. 1. How- 
ever, the result for the CHIK model is in better agreement 
with experiment with respect to the location of the min- 
imum and the density. As can be inferred from Fig. 1, 
at p = 2.2 g/cm 3 , a minimum value of about 0.35 GPa is 
found for the CHIK model. In contrast to this the BKS 
model has, at a density of p — 2.2 g/cm 3 , negative values 
for the pressure if T < 7000 K [19]. 

Figure 2 shows the partial pair correlation functions 
at T = 3600 K, as calculated from CPMD and classical 
MD using the BKS and the CHIK model. The CHIK 
model yields good agreement with the CPMD results. The 
largest differences are found in <?siSi( r )j but also in this 
case, the CHIK potential leads to a better agreement with 
CPMD than the BKS potential (see the inset in Fig. 2). 

The partial pair correlation functions as obtained from 
the CPMD run were used to parametrize the CHIK po- 
tential. That these functions are reproduced well by the 
CHIK potential is thus not that surprising. A less ob- 
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Fig. 2: Partial pair correlation functions at T = 3600 K, cal- 
culated from CPMD and classical MD using the BKS and the 
CHIK potentials, as indicated. The inset shows an enlargement 
of the first peak in <7SiSi(r)- 

vious test of the CHIK model is provided by consider- 
ing the distribution functions P Q /3 7 of the bond angles 
(a/37 = OSiO, SiOSi) in comparison to CPMD. As can be 
inferred from Fig. 3, the average intra-tctrahcdral OSiO 
angle is around 6 = 106° at 3600 K, both for CPMD and 
classical MD. However, in the case of the BKS model the 
function PosiO is less broad and exhibits a significantly 
higher amplitude than the CPMD and the CHIK model. 
Larger differences are seen in the distribution function for 
the inter-tetrahedral SiOSi angle. While at 3600 K the 
BKS distribution shows a maximum around 9 = 147°, the 
CHIK potential gives a maximum around 9 = 141°, in 
better agreement with the CPMD value around 9 = 136°. 

The distribution functions for the SiOSi angle exhibit 
also a shoulder around 9 = 90° revealing the emergence of 
edge-sharing tetrahedra at the relatively high temperature 
T = 3600 K [19]. In the BKS case, the latter shoulder has 
a lower amplitude, indicating a smaller number of edge- 
sharing tetrahedra. As for PoSiO, the BKS function for 
-PsiOSi is less broad and its main peak has a higher am- 
plitude. In this sense, the BKS model leads to a less dis- 
ordered structure than the CPMD and the CHIK model. 
Also included in Fig. 3 is the SiOSi bond angle distribu- 
tion at 300 K (i.e. in the glass), as obtained from the sim- 
ulation with the CHIK model (for the preparation of the 
glass samples, see below the discussion on the calculation 
of the density of states). In agreement with experimen- 
tal data [36-38] it shows a shift of the average value for 
the SiOSi angle from about 140° for the melt at 3600 K 
to about 150° for the glass structure (note the very good 
agreement with the value of 151°, obtained by an analysis 
of NMR data by Mauri et al [37]). In contrast to that 
the BKS potential shows only a slight change of the SiOSi 
angle from the melt to the glass structure [19]. 

Dynamic Features. The self-diffusion constants D a 
(a = Si, O) were computed from the long-time limit of 



0.040 J 




Fig. 3: Angular distribution functions for the OSiO and SiOSi 
angle at T = 3600 K. For the CHIK model, the SiOSi distribu- 
tion at T — 300 K is also included. The vertical lines indicate 
experimental values from Refs. [36] (dotted line), [37] (solid 
line), and [38] (dashed line). 

the mean squared displacements (r^(t)) via the Einstein 
relation D a = lim 4 ^ 00 (r^(t))/6i [31,32]. In Fig. 4, the 
temperature dependence of the D a is displayed in an Ar- 
rhenius plot. The diffusion dynamics as predicted by 
the CHIK model appears to be faster than that of the 
BKS model. At low temperatures (T ~ 2750 K), the self- 
diffusion coefficients are about a factor of 5 higher than 
those obtained from the BKS model. 

In agreement with various experimental studies 
[24], at low temperatures the self-diffusion constants 
can be well-described by an Arrhcnius law, D a = 
A a exp[—E"/(kBT)]. The activation energies E a , that 
we find for the CHIK model, are 4.51 eV for oxygen and 
4.97 eV for silicon. These values are very similar to those 
obtained for the BKS model [20] and they are in good 
agreement with experimental results (P a = 4.7 eV for oxy- 
gen [39] and P a = 6.0 eV for silicon [40]). 

To compute the vibrational density of states g{v) 
(with v the frequency), 16 independent and fully equili- 
brated samples at 2440 K were quenched instantaneously 
to 300 K, followed by an annealing for 80 ps. Then, g{u) 
was determined by calculating the Fourier transform of the 
velocity autocorrelation functions for Si and O (for details 
see Ref. [21]). In Fig. 5, the so-obtained density of states 
is compared to results from CPMD (from Ref. [10]) and 
MD simulations using the BKS potential (from Ref. [21]). 

A prominent feature in g{v) is the double-peak in the 
frequency band above about 28THz. The vibrational ex- 
citations in this band correspond to stretching modes of 
the Si-0 bond [31]. Note that the BKS potential has been 
fitted to reproduce the high frequency band of the vibra- 
tional spectrum. Thus, it is not surprising that it is in bet- 
ter agreement with CPMD than the CHIK model, since we 
have not included any vibrational properties in the fitting 
procedure of the CHIK potential. The CHIK model seems 
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Fig. 4: Arrhenius plot of the self-diffusion constants obtained 
by simulations with the BKS and the CHIK potentials. The 
bold solid lines are fits with Arrhenius laws (see text). 

to be better than the BKS model in the intermediate fre- 
quency band 20THz< v < 30THz where, in contrast to 
the BKS case, a single peak is observed, albeit at a slightly 
lower frequency than in CPMD result. Below 20THz the 
BKS and CHIK results are very similar and do not agree 
well with the CPMD result. In order to significantly im- 
prove the description of the density of states it might be 
necessary to account for polarization effects in the model 
potential, as suggested by Wilson et al. [16]. 

a— quartz. We have also carried out calculations on 
a— quartz to check the transferability of our potential. In 
the experimental phase diagram [31], a— quartz is the sta- 
ble phase in the low temperature domain (i.e. up to 846 K). 
As already mentioned, some of the crystalline properties of 
a— quartz, i.e. different lattice parameters and elastic con- 
stants, have been used for the parametrization of the BKS 
potential [2] . We have also calculated these quantities us- 
ing our potential (now using a cutoff at 10 A for the short- 
ranged part of the potential). For this, the GULP code [41] 
was employed which yields the different lattice properties 
at OK. The results are summarized in Tables (2) and (3), 
in comparison to those predicted by the BKS model and 
to experimental data (the latter at 300 K!). The CHIK po- 
tential reproduces the experimental data with comparable 
accuracy as the BKS potential, thus indicating that it is 
also well suited for the study of crystalline SiC>2 phases. 

Conclusions. — This work was devoted to the devel- 
opment of a new effective potential for silica. We have 
presented a fitting scheme for deriving effective potentials 
from ab initio simulations (CPMD). The new potential 
(named CHIK potential) is superior to the BKS model 
with respect to various static and dynamic properties of 
amorphous silica. The CHIK potential is also found to be 
transferable to other phases such as a-quartz, considering 
the quantitative agreement between the experimental data 
for both cell parameters and elastic constants. The fitting 




v (THz) 

Fig. 5: Density of states from the different MD simulations, as 
indicated. 

scheme proposed in this study is well-suited to parametrize 
potentials for other (amorphous) systems, in particular for 
mixtures of Si02 with other oxides such as alkali oxides or 
AI2O3 where the local structure of the melt or the glass 
can be very different from that of crystalline phases. 
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